/*************************************************************
*
* CAM_QUAL.CPP
*
* BVolume qualification against POV
*************************************************************/
#include <m3d/m3d_stdafx.h>

#include <ap_comm.h>
#include <m3d.h>
#include <camera.h>

//
// Static functions
//
//static int _camSeparateBBoxPOV(camPOV_QUAL *povQual,
//                               m3dV *axis, 
//                               m3dV *orgBB, m3dV axesBB[3],
//                               m3dV *orgCam, m3dV axesCam[4],
//                               float dCam);
static int _camSeparateOBBPOV (const camPOV_QUAL *povQual,
                               const m3dV *axis, 
                               const m3dV *orgBB, const m3dV axesBB[3],
                               const m3dV *orgCam, const m3dV axesCam[4],
                               float dCam);
//static int _camSeparateBBoxPOV_AxisPOV (camPOV_QUAL *povQual,
//                                        int axisNmb,
//                                        m3dV *orgBB, m3dV axesBB[3]);

/*************************************************************
* camCAMERA::CalcPOVQual()
*
* Compute data to perform speedy qualification by POV
* This function is to be called ater any changes to camera
* parameters, before rendering
*************************************************************/
void camCAMERA::CalcPOVQual (void)
{
   m3dV  v;
   int   i;

   // compute dCam
   this->povQual.dCam = m3dSqrt(m3dPow2(dxViewPort*0.5f) + 
                                m3dPow2(dyViewPort*0.5f) + 
                                m3dPow2(nearPlane));

   GetOrigin(&povQual.orgCam);

   // compute camera frustum points (with camera origin in (0;0;0))
   // OPTIMIZE: this can be done without matrix opertaions using up-left-forw vectors
   m3dV tmpPoints[8];
   for (i = 0; i < 4; i++)  {
      // near-plane points
      v.x = polyViewport.vert[i].x;
      v.y = polyViewport.vert[i].y;
      v.z = nearPlane;
      matrC2W.TransformVector(&v, &tmpPoints[i]);
      // far-plane points
      v.x *= farPlane / nearPlane;
      v.y *= farPlane / nearPlane;
      v.z  = farPlane;
      matrC2W.TransformVector(&v, &tmpPoints[i + 4]);
   }
   // compute frustum corner vectors
   for (i = 0; i < 4; i++)  {
      povQual.axisCam[i] = tmpPoints[i];
      m3dNormalize(&povQual.axisCam[i]);
   }
   // move camera frustum points to world-space
   for (i = 0; i < 8; i++)  {
      tmpPoints[i] += povQual.orgCam;
   }
   frustum.UpdatePoints(tmpPoints);
   // calculate normals to 4 side planes
   for (i = 0; i < 4; i++)  {
      m3dCrossProduct(&povQual.axisCam[i],
                      &povQual.axisCam[(i + 1) % 4],
                      &povQual.axisCamPerp[i]);
   }

   // compute the last camera-POV vector
   m3dMakeVZ(&v, -1.f);
   matrC2W.TransformVector(&v, &povQual.axisCamPerp[4]);

   povQual.dCamProd[0] = 0.f;
   povQual.dCamProd[1] = 0.f;
   povQual.dCamProd[2] = 0.f;
   povQual.dCamProd[3] = 0.f;

   // this is the constant to speed up qualification
   // by the 5th axis of POV
   povQual.dCamProd[4] = povQual.dCam*m3dDotProduct(&povQual.axisCam[0], &povQual.axisCamPerp[4]);

   // Note: reordered planes to correspond to F/L/R first, then bottom, then top, to optimize culling performance
   // calc frustum
   frustum.ConstructPlane(4, -povQual.axisCamPerp[0], povQual.orgCam);
   frustum.ConstructPlane(1, -povQual.axisCamPerp[1], povQual.orgCam);
   frustum.ConstructPlane(3, -povQual.axisCamPerp[2], povQual.orgCam);
   frustum.ConstructPlane(2, -povQual.axisCamPerp[3], povQual.orgCam);

   m3dV  camForw;
   GetVForw(&camForw);

   frustum.ConstructPlane(0, camForw, povQual.orgCam + camForw * nearPlane); // !!! Nearest plane must be 0 !!!
   frustum.SetNumPlanes(5);

   return;
}

/*************************************************************
* _camSeparateBBoxPOV()
*
*************************************************************/
//static int _camSeparateBBoxPOV (camPOV_QUAL *povQual,
//                                m3dV *axis, 
//                                m3dV *orgBB, m3dV axesBB[3],
//                                m3dV *orgCam, m3dV axesCam[4],
//                                float dCam)
//{
//   float prod[4], p;
//   int   sign;
//   float limPOV;
//   float centerBB, radBB;
//   float limBBUp, limBBDown;
//   //BOOL  signCalced = FALSE;
//
//   // ASSERT(!m3dIsZero(m3dLengthVector(axis)));
//   if (m3dIsZero_A(m3dLengthVector(axis), 0.001f))  {
//      return M3D_INSIDE;
//   }
//
//   prod[0] = m3dDotProduct(&axesCam[0],axis);
//   sign    = (prod[0] > M3D_EPSILON);
//
//   prod[1] = m3dDotProduct(&axesCam[1],axis);
//   if (sign ^ (prod[1] > M3D_EPSILON))  {
//      return M3D_INSIDE;
//   }
//
//   prod[2] = m3dDotProduct(&axesCam[2],axis);
//   if (sign ^ (prod[2] > M3D_EPSILON))  {
//      return M3D_INSIDE;
//   }
//
//   prod[3] = m3dDotProduct(&axesCam[3],axis);
//   if (sign ^ (prod[3] > M3D_EPSILON))  {
//      return M3D_INSIDE;
//   }
//
//   // now POV projection is limited from one side
//
//   if (sign > 0)  {
//      // POV projection is down-limited
//      p = m3dDotProduct(orgCam,axis);
//      limPOV = p + dCam*min(prod[0],min(prod[1],min(prod[2],prod[3])));
//      centerBB = m3dDotProduct(axis,orgBB);
//      radBB = m3dFAbs(axesBB[0].x*axis->x) + 
//              m3dFAbs(axesBB[1].y*axis->y) + 
//              m3dFAbs(axesBB[2].z*axis->z);
//      limBBUp = centerBB + radBB;
//      if (limPOV > limBBUp)  {
//         return M3D_OUTSIDE; // bbox is separated from POV along 'axis'
//      }
//      limBBDown = centerBB - radBB;
//      if (limBBDown > limPOV)  {
//         return M3D_INSIDE; // bbox is inside POV along 'axis'
//      }
//      return M3D_ISECT;    // bbox isects POV along 'axis'
//   }
//
//   ASSERT(sign == 0);
//
//   // POV projection is up-limited
//   p = m3dDotProduct(orgCam,axis);
//   limPOV = p + dCam*max(prod[0],max(prod[1],max(prod[2],prod[3])));
//   centerBB = m3dDotProduct(axis,orgBB);
//   radBB = m3dFAbs(axesBB[0].x*axis->x) + 
//           m3dFAbs(axesBB[1].y*axis->y) + 
//           m3dFAbs(axesBB[2].z*axis->z);
//   limBBDown = centerBB - radBB;
//   if (limBBDown > limPOV)  {
//      return M3D_OUTSIDE; // bbox is separated from POV along 'axis'
//   }
//   limBBUp = centerBB + radBB;
//   if (limPOV > limBBUp)  {
//      return M3D_INSIDE; // bbox is inside POV along 'axis'
//   }
//   return M3D_ISECT;     // bbox isects POV along 'axis'
//}

/*************************************************************
* _camSeparateOBBPOV()
*
*************************************************************/
static int _camSeparateOBBPOV (const camPOV_QUAL *povQual,
                               const m3dV *axis, 
                               const m3dV *orgBB, const m3dV axesBB[3],
                               const m3dV *orgCam, const m3dV axesCam[4],
                               float dCam)
{
   float prod[4], p;
   int   sign;
   float limPOV;
   float centerOBB, radOBB;
   float limOBBUp, limOBBDown;
   //BOOL  signCalced = FALSE;

   if (m3dIsZero_A(m3dLengthVector(axis), 0.001f))  {
      return M3D_INSIDE;
   }

   prod[0] = m3dDotProduct(&axesCam[0],axis);
   sign    = (prod[0] > M3D_EPSILON);

   prod[1] = m3dDotProduct(&axesCam[1],axis);
   if (sign ^ (prod[1] > M3D_EPSILON))  {
      return M3D_INSIDE;
   }

   prod[2] = m3dDotProduct(&axesCam[2],axis);
   if (sign ^ (prod[2] > M3D_EPSILON))  {
      return M3D_INSIDE;
   }

   prod[3] = m3dDotProduct(&axesCam[3],axis);
   if (sign ^ (prod[3] > M3D_EPSILON))  {
      return M3D_INSIDE;
   }

   // now POV projection is limited from one side

   if (sign > 0)  {
      // POV projection is down-limited
      p = m3dDotProduct(orgCam,axis);
      limPOV = p + dCam*min(prod[0],min(prod[1],min(prod[2],prod[3])));
      centerOBB = m3dDotProduct(axis,orgBB);
      radOBB = m3dFAbs(m3dDotProduct(&axesBB[0],axis)) + 
               m3dFAbs(m3dDotProduct(&axesBB[1],axis)) + 
               m3dFAbs(m3dDotProduct(&axesBB[2],axis));
      limOBBUp = centerOBB + radOBB;
      if (limPOV > limOBBUp)  {
         return M3D_OUTSIDE; // bbox is separated from POV along 'axis'
      }
      limOBBDown = centerOBB - radOBB;
      if (limOBBDown > limPOV)  {
         return M3D_INSIDE; // bbox is inside POV along 'axis'
      }
      return M3D_ISECT;    // bbox isects POV along 'axis'
   }

   ASSERT(sign == 0);

   // POV projection is up-limited
   p = m3dDotProduct(orgCam,axis);
   limPOV = p + dCam*max(prod[0],max(prod[1],max(prod[2],prod[3])));
   centerOBB = m3dDotProduct(axis,orgBB);
   radOBB = m3dFAbs(m3dDotProduct(&axesBB[0],axis)) + 
            m3dFAbs(m3dDotProduct(&axesBB[1],axis)) + 
            m3dFAbs(m3dDotProduct(&axesBB[2],axis));
   limOBBDown = centerOBB - radOBB;
   if (limOBBDown > limPOV)  {
      return M3D_OUTSIDE; // bbox is separated from POV along 'axis'
   }
   limOBBUp = centerOBB + radOBB;
   if (limPOV > limOBBUp)  {
      return M3D_INSIDE; // bbox is inside POV along 'axis'
   }
   return M3D_ISECT;     // bbox isects POV along 'axis'
}

/*************************************************************
* _camSeparateBBoxPOV_AxisPOV()
*
*************************************************************/
//static int _camSeparateBBoxPOV_AxisPOV (camPOV_QUAL *povQual,
//                                        int axisNmb,
//                                        m3dV *orgBB, m3dV axesBB[3])
//{
//   float limPOV;
//   float centerBB, radBB;
//   float limBBUp, limBBDown;
//   m3dV  axis;
//
//   axis = povQual->axisCamPerp[axisNmb];
//
//   // POV projection is always zero-up-limited
//   limPOV = m3dDotProduct(&povQual->orgCam,&axis) + povQual->dCamProd[axisNmb];
//   centerBB = m3dDotProduct(&axis,orgBB);
//   radBB = m3dFAbs(axesBB[0].x*axis.x) + 
//           m3dFAbs(axesBB[1].y*axis.y) + 
//           m3dFAbs(axesBB[2].z*axis.z);
//   limBBDown = centerBB - radBB;
//   if (limBBDown > limPOV)  {
//      return M3D_OUTSIDE; // bbox is separated from POV along 'axis'
//   }
//   limBBUp = centerBB + radBB;
//   if (limPOV > limBBUp)  {
//      return M3D_INSIDE; // bbox is inside POV along 'axis'
//   }
//   return M3D_ISECT;     // bbox isects POV along 'axis'
//}

/*************************************************************
* _camSeparateOBBPOV_AxisPOV()
*
*************************************************************/
static int _camSeparateOBBPOV_AxisPOV (const camPOV_QUAL *povQual,
                                       int axisNmb, 
                                       const m3dV *orgBB, m3dV axesBB[3])
{
   float limPOV;
   float centerOBB, radOBB;
   float limOBBUp, limOBBDown;
   m3dV  axis;

   axis = povQual->axisCamPerp[axisNmb];

   if (m3dIsZero_A(m3dLengthVector(&axis), 0.001f))  {
      return M3D_INSIDE;
   }

   // POV projection is always zero-up-limited
   limPOV = m3dDotProduct(&povQual->orgCam,&axis) + povQual->dCamProd[axisNmb];
   centerOBB = m3dDotProduct(&axis,orgBB);
   radOBB = m3dFAbs(m3dDotProduct(&axesBB[0],&axis)) + 
            m3dFAbs(m3dDotProduct(&axesBB[1],&axis)) + 
            m3dFAbs(m3dDotProduct(&axesBB[2],&axis));
   limOBBDown = centerOBB - radOBB;
   if (limOBBDown > limPOV)  {
      return M3D_OUTSIDE; // bbox is separated from POV along 'axis'
   }
   limOBBUp = centerOBB + radOBB;
   if (limPOV > limOBBUp)  {
      return M3D_INSIDE; // bbox is inside POV along 'axis'
   }
   return M3D_ISECT;     // bbox isects POV along 'axis'
}




/*************************************************************
* camCAMERA::QualOBBPOV()
* 
*************************************************************/
int camCAMERA::QualOBBPOV (const m3dOBB *pOBB) const
{
   m3dV  axisOBB[3];    // 3 obb half-axes
   m3dV  orgOBB;        // obb origin
   m3dV  axis;
   int   i, j;
   int   nInside;       // number of 'inside' qualifications
   int   rc;

   // compute obb axes/origin (recompute only non-zero components)
   m3dScaleVector(pOBB->GetDirX(), 0.5f * pOBB->GetSizeX(), &axisOBB[0]);
   m3dScaleVector(pOBB->GetDirY(), 0.5f * pOBB->GetSizeY(), &axisOBB[1]);
   m3dScaleVector(pOBB->GetDirZ(), 0.5f * pOBB->GetSizeZ(), &axisOBB[2]);

   pOBB->GetCenter(&orgOBB);

   nInside = 0; // count the number of 'inside' qualifications

   // separate by POV axes
   for (i = 0; i < 5; i++)  {
      rc = _camSeparateOBBPOV_AxisPOV(&this->povQual,i, &orgOBB, axisOBB);
      ASSERT(rc == _camSeparateOBBPOV(&this->povQual,
      &this->povQual.axisCamPerp[i], 
      &orgOBB, axisOBB,
      &this->povQual.orgCam, this->povQual.axisCam,
      this->povQual.dCam));
      if (rc == M3D_OUTSIDE)  {
         return M3D_OUTSIDE;
      }
      if (rc == M3D_INSIDE)  {
         nInside++;
      }
   }

   if (nInside == 5)  {
      // this is the same as the qualification
      // by 5 coordinate planes of POV
      return M3D_INSIDE;
   }

   // separate by OBB axes
   for (i = 0; i < 3; i++)  {
      rc = _camSeparateOBBPOV(&this->povQual,
            &axisOBB[i], 
            &orgOBB, axisOBB,
            &this->povQual.orgCam, this->povQual.axisCam,
            this->povQual.dCam);
      if (rc == M3D_OUTSIDE)  {
         return M3D_OUTSIDE;
      }
      if (rc == M3D_INSIDE)  {
         nInside++;
      }
   }

   // separate by cross-product axes
   for (i = 0; i < 3; i++)  {
      for (j = 0; j < 5; j++)  {
         m3dCrossProduct(&axisOBB[i], &this->povQual.axisCamPerp[j], &axis);
         rc = _camSeparateOBBPOV(&this->povQual,
               &axis, 
               &orgOBB, axisOBB,
               &this->povQual.orgCam, this->povQual.axisCam,
               this->povQual.dCam);
         if (rc == M3D_OUTSIDE)  {
            return M3D_OUTSIDE;
         }
         if (rc == M3D_INSIDE)  {
            nInside++;
         }
      }
   }

   if (nInside == 23)  {
      return M3D_INSIDE;
   }

   return M3D_ISECT;
}
/*************************************************************
* camCAMERA::QualBoxPOV ()
* 
*************************************************************/
int camCAMERA::QualBoxPOV (const m3dBOX *pBB)
{
   m3dOBB obb;
   pBB->ConvToOBB(&obb);
   return QualOBBPOV(&obb);
}


/*************************************************************
* camCAMERA::QualPointPOV()
*
*************************************************************/
int camCAMERA::QualPointPOV (const m3dV *p)
{
   m3dV     vCCS;
   float    ooz;
   m2dPOLY *poly;

   TransformPoint(p, &vCCS);
   if (vCCS.z <= this->nearPlane) {
      return M3D_OUTSIDE;
   }
   // project to image plane
   ooz = nearPlane / vCCS.z;
   vCCS.x *= ooz;
   vCCS.y *= ooz;

   poly = &this->polyViewport;
   if (vCCS.x < m3dPolyVertexX(poly, 0) || vCCS.y < m3dPolyVertexY(poly, 0) || 
       vCCS.x > m3dPolyVertexX(poly, 1) || vCCS.y > m3dPolyVertexY(poly, 3)) {
      return M3D_OUTSIDE;
   }
   return M3D_INSIDE;
}
/*************************************************************
* camCAMERA::QualSpherePOV()
* 
* Algorithm computes closest point in all cases by determining in which
* Voronoi region of the vertices, edges, and faces of the frustum that
* the test point lives.
*************************************************************/
int camCAMERA::QualSpherePOV (const m3dSPHERE *pPS)
{
   m3dV  vD, vT, vC;  
   m3dV  center;
   
   pPS->GetCenter(&center);
   m3dSubtractVector(&center, &povQual.orgCam, &vD);
   
   // Transform point to frustum CS
   vT.x = m3dDotProduct(&vD, GetVRight());
   vT.y = m3dDotProduct(&vD, GetVUp());
   vT.z = m3dDotProduct(&vD, GetVForw());

   // move point to positive side 
   if (vT.x < 0.f) vT.x = -vT.x;
   if (vT.y < 0.f) vT.y = -vT.y;

   float depthRatio = farPlane / nearPlane; 
   float lMin = 0.5f * GetDXViewPort();
   float lMax = depthRatio * lMin;
   float uMin = 0.5f * GetDYViewPort();
   float uMam = depthRatio * uMin;
   float dMin = nearPlane;
   float dMax = farPlane;
   float lMin_2 = m3dPow2(lMin);
   float uMin_2 = m3dPow2(uMin);
   float dMin_2 = m3dPow2(dMin);
   float ldDotMin  = lMin_2 + dMin_2;
   float udDotMin  = uMin_2 + dMin_2;
   float ludDotMin = lMin_2 + udDotMin;
   float ldDotMax  = depthRatio * ldDotMin;
   float udDotMax  = depthRatio * udDotMin;
   float ludDotMax = depthRatio * ludDotMin;

   float lDot, uDot, ldDot, udDot, ludDot, lEdgeDot, uEdgeDot, t;

   if ( vT.z >= dMax ) {
      if ( vT.x <= lMax ) {
         if (vT.y <= uMam) { vC.x = vT.x; vC.y = vT.y; vC.z = dMax; } 
            else { vC.x = vT.x; vC.y = uMam; vC.z = dMax; }
      } else   {
         if (vT.y <= uMam) { vC.x = lMax;   vC.y = vT.y; vC.z = dMax; } 
            else { vC.x = lMax;   vC.y = uMam; vC.z = dMax; }
      }
   } else if (vT.z <= dMin)  {
      if (vT.x <= lMin)  {
         if (vT.y <= uMin) { 
            vC.x = vT.x; vC.y = vT.y; vC.z = dMin; 
         } else {
            udDot = uMin * vT.y + dMin * vT.z;
            if (udDot >= udDotMax) { vC.x = vT.x; vC.y = uMam; vC.z = dMax; } 
               else if (udDot >= udDotMin) { uDot = dMin * vT.y - uMin * vT.z; t = uDot/udDotMin; vC.x = vT.x; vC.y = vT.y - t * dMin; vC.z = vT.z + t * uMin;} 
               else { vC.x = vT.x; vC.y = uMin; vC.z = dMin; }
         }
      } else   {
         if (vT.y <= uMin)  {
            ldDot = lMin * vT.x + dMin * vT.z;
            if (ldDot >= ldDotMax) { vC.x = lMax; vC.y = vT.y; vC.z = dMax; } 
               else if (ldDot >= ldDotMin) { lDot = dMin * vT.x - lMin * vT.z; t = lDot/ldDotMin; vC.x = vT.x - t * dMin; vC.y = vT.y; vC.z = vT.z + t * lMin; } 
               else { vC.x = lMin; vC.y = vT.y; vC.z = dMin; }
         } else   {
            ludDot = lMin*vT.x + uMin*vT.y + dMin*vT.z;
            lEdgeDot = uMin * ludDot - ludDotMin * vT.y;
            if ( lEdgeDot >= 0.f ) {
               ldDot = lMin*vT.x + dMin*vT.z;
               if (ldDot >= ldDotMax) { vC.x = lMax; vC.y = vT.y; vC.z = dMax; } 
                  else if (ldDot >= ldDotMin) { lDot = dMin*vT.x - lMin*vT.z; t = lDot/ldDotMin; vC.x = vT.x - t * dMin; vC.y = vT.y; vC.z = vT.z + t * lMin; } 
                  else { vC.x = lMin; vC.y = vT.y; vC.z = dMin; }
            } else   {
               uEdgeDot = lMin * ludDot - ludDotMin * vT.x;
               if ( uEdgeDot >= 0.f ) {
                  udDot = uMin*vT.y + dMin*vT.z;
                  if (udDot >= udDotMax) { vC.x = vT.x; vC.y = uMam; vC.z = dMax; } 
                     else if (udDot >= udDotMin) { uDot = dMin*vT.y - uMin*vT.z; t = uDot/udDotMin; vC.x = vT.x; vC.y = vT.y - t * dMin; vC.z = vT.z + t * uMin; } 
                     else { vC.x = vT.x; vC.y = uMin; vC.z = dMin; }
               } else   {
                  if (ludDot >= ludDotMax) { vC.x = lMax; vC.y = uMam; vC.z = dMax; } 
                     else if (ludDot >= ludDotMin) { t = ludDot/ludDotMin; vC.x = t * lMin; vC.y = t * uMin; vC.z = t * dMin; } 
                     else { vC.x = lMin; vC.y = uMin; vC.z = dMin; }
               }
            }
         }
      }
   } else   {
      lDot = dMin * vT.x - lMin * vT.z;
      uDot = dMin * vT.y - uMin * vT.z;
      if (lDot <= 0.f)  {
         if (uDot <= 0.f)  {
            vC = vT;
         } else   {
            udDot = uMin * vT.y + dMin * vT.z;
            if (udDot >= udDotMax) { vC.x = vT.x; vC.y = uMam; vC.z = dMax; } 
               else { uDot = dMin * vT.y - uMin * vT.z; t = uDot/udDotMin; vC.x = vT.x; vC.y = vT.y - t * dMin; vC.z = vT.z + t * uMin; }
         }
      } else   {
         if (uDot <= 0.f)  {
            ldDot = lMin * vT.x + dMin * vT.z;
            if (ldDot >= ldDotMax) { vC.x = lMax; vC.y = vT.y; vC.z = dMax; } 
               else { lDot = dMin * vT.x - lMin * vT.z; t = lDot/ldDotMin; vC.x = vT.x - t * dMin; vC.y = vT.y; vC.z = vT.z + t * lMin; }
         } else   {
            ludDot = lMin * vT.x + uMin * vT.y + dMin * vT.z;
            lEdgeDot = uMin * ludDot - ludDotMin * vT.y;
            if (lEdgeDot >= 0.f) {
               ldDot = lMin * vT.x + dMin * vT.z;
               if (ldDot >= ldDotMax) { vC.x = lMax; vC.y = vT.y; vC.z = dMax; } 
                  else { lDot = dMin * vT.x - lMin * vT.z; t = lDot/ldDotMin; vC.x = vT.x - t * dMin; vC.y = vT.y; vC.z = vT.z + t * lMin; }
            } else   {
               uEdgeDot = lMin * ludDot - ludDotMin * vT.x;
               if (uEdgeDot >= 0.f) {
                  udDot = uMin * vT.y + dMin * vT.z;
                  if (udDot >= udDotMax) { vC.x = vT.x; vC.y = uMam; vC.z = dMax; } 
                     else { uDot = dMin * vT.y - uMin * vT.z; t = uDot/udDotMin; vC.x = vT.x; vC.y = vT.y - t * dMin; vC.z = vT.z + t * uMin; }
               } else   {
                  if (ludDot >= ludDotMax)  { vC.x = lMax; vC.y = uMam; vC.z = dMax; } 
                     else { t = ludDot/ludDotMin; vC.x = t * lMin; vC.y = t * uMin; vC.z = t * dMin; }
               }
            }
         }
      }
   }
   m3dSubtractVector(&vT, &vC, &vD);

   if (m3dLengthVector_2(&vD) <= pPS->GetRadius() * pPS->GetRadius()) {
      return M3D_INSIDE;
   }
   return M3D_OUTSIDE;
}

/*************************************************************
* m3dQualSpherePOV()
*
* IN: org     - sphere origin
*     radius  - radius of the sphere
* 
* Qualify BSphere by POV. This function intersects cone 
* (sphere is seen by a cone) with plane of projection, 
* builds bbox of ellipse and qualifies it by viewport.
*************************************************************/
/*
int camCAMERA::QualSpherePOV (m3dSPHERE *pPS)
{
   m3dV     orgCCS;
   m2dPOLY  bbox2D, bboxClip2D;
   float    cos2_alpha;
   float    x1, x2;
   float    y1, y2;
   float    a, b, c, d;
   float    coeff1, coeff2, coeff3;
   int      nroot;
   m3dV    *org = &pPS->org;
   float    radius = pPS->radius;

   ASSERT(radius  > M3D_EPSILON);
   TransformPoint(org, &orgCCS);

   if (orgCCS.z + radius < this->nearPlane)  {
      return M3D_OUTSIDE;
   } 
   float dist = m3dLengthVector(&orgCCS);
   if (dist - radius < this->nearPlane+100.f*M3D_EPSILON)  {   //orgCCS.z  - error
      return M3D_ISECT;
   }

   cos2_alpha = 1.f - m3dPow2(radius/m3dLengthVector(&orgCCS));
   // check the range of cos2_alpha
   ASSERT(cos2_alpha > -M3D_EPSILON);
   ASSERT(cos2_alpha < 1.f  + M3D_EPSILON);

   m3dNormalize(&orgCCS);
   
   a = orgCCS.x;
   b = orgCCS.y;
   c = orgCCS.z;
   d = nearPlane;

   coeff1 = a*a+b*b-cos2_alpha;
   coeff2 = 2*c*d;
   coeff3 = c*c-cos2_alpha;
   nroot = m3dSolveQuadraticEq(coeff1,coeff2*b,d*d*(a*a+coeff3),&y1,&y2);
   ASSERT(nroot > 0);
   ASSERT(y1 >= y2);     // because a*a+b*b-cos2_alpha < 0

   nroot = m3dSolveQuadraticEq(coeff1,coeff2*a,d*d*(b*b+coeff3),&x1,&x2);
   ASSERT(nroot > 0);
   ASSERT(x1 >= x2);     // because a*a+b*b-cos2_alpha < 0

   m3dPolyNVert(&bbox2D) = 4;
   m3dPolyVertexX(&bbox2D,0) = x2;
   m3dPolyVertexY(&bbox2D,0) = y1;
   m3dPolyVertexX(&bbox2D,1) = x2;
   m3dPolyVertexY(&bbox2D,1) = y2;
   m3dPolyVertexX(&bbox2D,2) = x1;
   m3dPolyVertexY(&bbox2D,2) = y2;
   m3dPolyVertexX(&bbox2D,3) = x1;
   m3dPolyVertexY(&bbox2D,3) = y1;

   return m2dClipPolyRect(&bbox2D, &this->polyViewport, &bboxClip2D);
}
*/
/*************************************************************
* camCAMERA::QualBoxPOV()
*
* Qualification based on separating axis theorem - fastest
* release
*************************************************************/
/*
int camCAMERA::QualBoxPOV (m3dBOX *pBB)
{
   m3dV  axisBB[3] = {m3dVZero, m3dVZero, m3dVZero}; // 3 bbox half-axes
   m3dV  orgBB;            // bbox origin
   m3dV  axis;
   int   i, j;
   int   nInside;          // number of 'inside' qualifications
   int   rc;

   // compute bbox axes/origin (recompute only non-zero components)
   axisBB[0].x = (pBB->fur.x - pBB->bll.x) * 0.5f;
   axisBB[1].y = (pBB->fur.y - pBB->bll.y) * 0.5f;
   axisBB[2].z = (pBB->fur.z - pBB->bll.z) * 0.5f;

   orgBB.x = (pBB->fur.x + pBB->bll.x) * 0.5f;
   orgBB.y = (pBB->fur.y + pBB->bll.y) * 0.5f;
   orgBB.z = (pBB->fur.z + pBB->bll.z) * 0.5f;

   // perform tests w/o optimizations
   nInside = 0; // count the number of 'inside' qualifications

   // separate by POV axes
   for (i = 0; i < 5; i++)  {
      rc = _camSeparateBBoxPOV_AxisPOV(&this->povQual,i,&orgBB,axisBB);
      if (rc == M3D_OUTSIDE)  {
         return M3D_OUTSIDE;
      }
      if (rc == M3D_INSIDE)  {
      nInside++;
      }

      ASSERT(rc == _camSeparateBBoxPOV(&this->povQual,
      &this->povQual.axisCamPerp[i], 
      &orgBB,axisBB,
      &this->povQual.orgCam,this->povQual.axisCam,
      this->povQual.dCam));
   }

   if (nInside == 5)  {
      // this is the same as the qualification
      // by 5 coordinate planes of POV
      return M3D_INSIDE;
   }

   // separate by BBox axes
   for (i = 0; i < 3; i++)  {
      rc = _camSeparateBBoxPOV(&this->povQual,
            &axisBB[i], 
            &orgBB,axisBB,
            &this->povQual.orgCam,this->povQual.axisCam,
            this->povQual.dCam);
      if (rc == M3D_OUTSIDE)  {
         return M3D_OUTSIDE;
      }
      if (rc == M3D_INSIDE)  {
         nInside++;
      }
   }

   // separate by cross-product axes
   for (i = 0; i < 3; i++)  {
      for (j = 0; j < 5; j++)  {
         m3dCrossProduct(&axisBB[i],&this->povQual.axisCamPerp[j],&axis);
         rc = _camSeparateBBoxPOV(&this->povQual,
               &axis, 
               &orgBB,axisBB,
               &this->povQual.orgCam,this->povQual.axisCam,
               this->povQual.dCam);
         if (rc == M3D_OUTSIDE)  {
            return M3D_OUTSIDE;
         }
         if (rc == M3D_INSIDE)  {
            nInside++;
         }
      }
   }

   if (nInside == 23)  {
      return M3D_INSIDE;
   }

   return M3D_ISECT;
}
*/

//
// End-of-file CAM_QUAL.CPP
//
